
// #define imax
#define cv(i, j)  cv[(i)*3 + j]
#define p(i)  p[i]
#define area(i) area[i]


__kernel void 
InitFlow(global double* cv, global double* p,  
        global const double* area, const double rho, 
        const double massFlow, const double p2,
        const double E)
{
    size_t gid = get_global_id(0);
    cv(gid, 0) = rho * area(gid);
    cv(gid, 1) = massFlow;
    cv(gid, 2) = rho * E * area(gid);
    p(gid) = p2;
}
// 宏定义有风险   只是原样替换，小心
#define du(i, j) du[(i)*3 + j]
#define rs(i, j) rs[(i)*3 + j]
#define ls(i, j) ls[(i)*3 + j]
#define vol(i)   vol[i]


inline double MUSCL3(double af, double bf, double eps)
{
    return (bf * (2.0 * af * af + eps) + af * (bf * bf + 2.0 * eps)) /
            (2.0 * af * af + 2.0 * bf * bf - af * bf + 3.0 * eps + 1.0e-30);
}

kernel void 
LrStateRoe(global const double*cv, global const double*area,
            global double* du, global double* rs, 
            global double* ls, global const double* p,
            global const double *vol,
            const int iorder, const double eps2rho,
            const double eps2u, const double eps2p)
{
    // 全局工作项有 ib2 个
    size_t gid = get_global_id(0);
    // size_t gidp1 = gid+1;
    double deltr[3], deltl[3];

    if(2 <= iorder)
    {
        du(gid+1, 0) = cv(gid+1, 0) / area( gid+1) - cv(gid, 0) / area(gid);
        du(gid+1, 1) = cv(gid+1, 1) / cv(gid+1, 0) - cv(gid, 1) / cv(  gid, 0);
        du(gid+1, 2) = p( gid+1) - p(gid); 

        // du[gidp1*3]  = cv[gidp1*3]/area[gidp1] - cv[gid*3]/area[gid];
        // du[gidp1*3 + 1] = cv[gidp1*3 +1]/cv[gidp1*3] -cv[gid*3 + 1]/cv[gid*3];
        // du[gidp1*3 + 2] = p[gidp1] - p[gid];
        if(gid == 0)
        {
            du(0, 0) = du(1, 0);
            du(0, 1) = du(1, 1);
            du(0, 2) = du(1, 2);
        }
        // 是 ib2-1 不是 ib2 注意
        if(gid == (ib2-1))
        {
            du(imax, 0) = du(ib2, 0);
            du(imax, 1) = du(ib2, 1);
            du(imax, 2) = du(ib2, 2);
        }
    }
    barrier(CLK_LOCAL_MEM_FENCE);

    double volt = 0.0;
    if(iorder == 3)
    {
        volt = pow(0.5 * (vol(gid + 1) + vol(gid)), 1.5);
        deltr[0] = MUSCL3(du(gid + 2, 0), du(gid + 1, 0), eps2rho * volt);
        deltl[0] = MUSCL3(du(gid + 1, 0), du(gid, 0), eps2rho * volt);

        deltr[1] = MUSCL3(du(gid + 2, 1), du(gid + 1, 1), eps2u * volt);
        deltl[1] = MUSCL3(du(gid + 1, 1), du(gid, 1), eps2u * volt);

        deltr[2] = MUSCL3(du(gid + 2, 2), du(gid + 1, 2), eps2p * volt);
        deltl[2] = MUSCL3(du(gid + 1, 2), du(gid, 2), eps2p * volt);

        rs(gid, 0) = cv(gid + 1, 0) / area(gid + 1) - 0.5 * deltr[0];
        rs(gid, 1) = cv(gid + 1, 1) / cv(gid + 1, 0) - 0.5 * deltr[1];
        rs(gid, 2) = p(gid + 1) - 0.5 * deltr[2];

        ls(gid, 0) = cv(gid, 0) / area(gid) + 0.5 * deltl[0];
        ls(gid, 1) = cv(gid, 1) / cv(gid, 0) + 0.5 * deltl[1];
        ls(gid, 2) = p(gid) + 0.5 * deltl[2];
    }
    // barrier(CLK_LOCAL_MEM_FENCE);
}


#define flux(i, j) flux[(i)*3 + j]
#define rhs(i, j)  rhs[(i)*3 + j]

inline double EntropyCorrection(double lam, double dlt)
{
    return (dlt < lam) ? lam : (lam * lam + dlt * dlt) / dlt * 0.5;
}
inline double mabs(double a)
{
    return (0 < a) ? a : -a;
}

kernel void
FluxRoeKer(global const double* area, global const double* ls,
            global const double* rs, global double* flux,
            global double* rhs, const double gamma,
            const double epsentr)
{
    // ib2 个工作项
    size_t gid = get_global_id(0);
    double gam1 = gamma - 1.0;
    double avs = 0.0;
    double avflux[3], lrho, lu, lp, lq, lh, rrho, ru, rp, rq, rh;
    double rhot, ut, ht, ct, q2t, c2t, tmp, fdiss[3];
    double lambda[3], delta, dv, eps1, eps2, eps3;

    // 面积为单元中心的面积，需要求个平均
    avs = 0.5 * (area(gid) + area(gid + 1));

    lrho = ls(gid, 0);
    lu = ls(gid, 1);
    lp = ls(gid, 2);
    lq = lu * lrho;
    // 总颔
    lh = gamma / gam1 * lp / lrho + 0.5 * lu * lu;

    rrho = rs(gid, 0);
    ru = rs(gid, 1);
    rp = rs(gid, 2);
    rq = ru * rrho;
    rh = gamma / gam1 * rp / rrho + 0.5 * ru * ru;
    // average of left and right convective fluxes
    avflux[0] = lq + rq;
    avflux[1] = lq * lu + lp + rq * ru + rp;
    avflux[2] = lq * lh + rq * rh;

    // dissipative fluxes
    rhot = sqrt(lrho * rrho);
    tmp = rhot / lrho;
    ut = (lu + ru * tmp) / (1.0 + tmp);
    ht = (lh + rh * tmp) / (1.0 + tmp);
    q2t = ut * ut;
    c2t = gam1 * (ht - 0.5 * q2t);
    ct = sqrt(c2t);

    lambda[0] = mabs(ut - ct);
    lambda[1] = mabs(ut);
    lambda[2] = mabs(ut + ct);
    delta = epsentr * ct;
    lambda[0] = EntropyCorrection(lambda[0], delta);
    lambda[1] = EntropyCorrection(lambda[1], delta);
    lambda[2] = EntropyCorrection(lambda[2], delta);
    dv = ru - lu;
    eps1 = lambda[0] * (rp - lp - rhot * ct * dv) / c2t * 0.5;
    eps2 = lambda[1] * (rrho - lrho - (rp - lp) / c2t);
    eps3 = lambda[2] * (rp - lp + rhot * ct * dv) / c2t * 0.5;

    fdiss[0] = eps1 + eps2 + eps3;
    fdiss[1] = eps1 * (ut - ct) + eps2 * ut + eps3 * (ut + ct);
    fdiss[2] = eps1 * (ht - ct * ut) + eps2 * q2t * 0.5 + eps3 * (ht + ct * ut);
    // 面总的flux
    flux(gid, 0) = 0.5 * avs * (avflux[0] - fdiss[0]);
    flux(gid, 1) = 0.5 * avs * (avflux[1] - fdiss[1]);
    flux(gid, 2) = 0.5 * avs * (avflux[2] - fdiss[2]);

    barrier(CLK_LOCAL_MEM_FENCE);

    if(gid != 0)
    {
        rhs(gid, 0) = flux(gid, 0) - flux(gid - 1, 0);
        rhs(gid, 1) = flux(gid, 1) - flux(gid - 1, 1);
        rhs(gid, 2) = flux(gid, 2) - flux(gid - 1, 2);
    }

}